#!/usr/local/bin/python
# Author: Jimmy Saw
# Date: 12/29/2011
# Description: This script extracts conflict regions (~100bp)
# Usage: resolve_conflicts.py seq.fasta sanitycheck.file

import sys
import re
from Bio import SeqIO
from Bio.SeqRecord import SeqRecord

seqfile = sys.argv[1]
sequence = SeqIO.read(seqfile, "fasta")

scfile = open(sys.argv[2], "rU")
sf = scfile.readlines()

for line in sf:
    tmp = line.split('\t')
    rn = tmp[0].zfill(2)
    loc = int(tmp[1])
    toextract = sequence.seq[loc-50:loc+50]
    tmprec = SeqRecord(toextract, id="conflict_"+rn, 
        description=str(loc-50)+":"+str(loc+50))
    print tmprec.format("fasta").strip('\n')

scfile.close()
